1.0 Prerequesits

# ReadLoad libraries and load shapefile with sf package
library(tidyverse)
library(sf)
library(raster)
library(leaflet)
setwd("~/Github/dataStories/topoSpec")
zoning <- st_read("zoning_subset_wgs84.shp")
Reading layer `zoning_subset_wgs84' from data source `/home/ryangarnett/Github/dataStories/topoSpec/zoning_subset_wgs84.shp' using driver `ESRI Shapefile'
Simple feature collection with 201 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -79.63432 ymin: 43.5901 xmax: -79.11799 ymax: 43.81809
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
pdm <- st_read("pdm_subset_wgs84.shp")
Reading layer `pdm_subset_wgs84' from data source `/home/ryangarnett/Github/dataStories/topoSpec/pdm_subset_wgs84.shp' using driver `ESRI Shapefile'
Simple feature collection with 12 features and 20 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -79.38804 ymin: 43.73857 xmax: -79.34854 ymax: 43.77281
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
pdm <- select(pdm, geometry)
head(pdm)
Simple feature collection with 6 features and 0 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -79.38804 ymin: 43.73857 xmax: -79.35253 ymax: 43.77281
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                        geometry
1 POLYGON ((-79.38417 43.7480...
2 POLYGON ((-79.3823 43.73857...
3 POLYGON ((-79.36177 43.7608...
4 POLYGON ((-79.3619 43.76225...
5 POLYGON ((-79.38592 43.7568...
6 POLYGON ((-79.37211 43.7408...
# Intersect zoning with pdm boundaries
zoning_int <- st_intersection(zoning, pdm)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries
count(zoning_int$zn_lclass)
head(zoning_int)
Simple feature collection with 6 features and 5 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -79.3849 ymin: 43.74801 xmax: -79.38146 ymax: 43.75797
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
   OBJECTID ZN_ZONE zn_class                zn_desc zn_lclass                       geometry
8      9017      ON        O     Open Space Natural     Rural MULTIPOLYGON (((-79.38471 4...
11    15008      RD        R   Residential Detached  Suburban POLYGON ((-79.3849 43.75182...
12     8144      RD        R   Residential Detached  Suburban POLYGON ((-79.38485 43.7515...
41     6154      OR        O  Open Space Recreation     Rural MULTIPOLYGON (((-79.38414 4...
46     5160      RD        R   Residential Detached  Suburban POLYGON ((-79.38249 43.7578...
52    14814      CR        C Commercial Residential      Core POLYGON ((-79.38414 43.7480...
plot(st_geometry(zoning_int), col = sf.colors(4, categorical = TRUE), border = 'grey', axes = TRUE)

#plot(st_geometry(nc), col = sf.colors(12, categorical = TRUE), border = 'grey', 
#    axes = TRUE)
zoning_int <- mutate(zoning_int, zn_num = ifelse(zn_lclass == "Core", 1,
                      ifelse(zn_lclass == "Industrial", 2,
                      ifelse(zn_lclass == "Rural", 3,
                      ifelse(zn_lclass == "Suburban", 4,
                      0)))))
head(zoning_int)
Simple feature collection with 6 features and 26 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -79.3849 ymin: 43.74801 xmax: -79.38146 ymax: 43.75797
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
   OBJECTID ZN_ZONE zn_class                zn_desc zn_lclass OBJECTID.1  AREA_ID    PDM NEXTPLANNE LASTUPDATE UPDATEDBY HYPERLINK  TIMESTAMP UPDATEUSER        X       Y
8      9017      ON        O     Open Space Natural     Rural        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
11    15008      RD        R   Residential Detached  Suburban        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
12     8144      RD        R   Residential Detached  Suburban        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
41     6154      OR        O  Open Space Recreation     Rural        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
46     5160      RD        R   Residential Detached  Suburban        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
52    14814      CR        C Commercial Residential      Core        557 20136416 52N-11       <NA>       2005      <NA>      <NA> 2015-11-12    GCCGDIA 314381.8 4845727
   LONGITUDE LATITUDE SHAPE_AREA SHAPE_LEN X2016 area_km ecs_2016 PNTCNT area_raw zn_num                       geometry
8  -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      3 MULTIPOLYGON (((-79.38471 4...
11 -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      4 POLYGON ((-79.3849 43.75182...
12 -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      4 POLYGON ((-79.38485 43.7515...
41 -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      3 MULTIPOLYGON (((-79.38414 4...
46 -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      4 POLYGON ((-79.38249 43.7578...
52 -79.38082 43.75343  0.6673742         0     0    0.67     <NA>      4 667374.2      1 POLYGON ((-79.38414 43.7480...
palette <- colorBin(c('#7b3294',
                      '#c2a5cf',
                      '#a6dba0',
                      '#008837'), 
                     bins = c(1, 2, 3, 4))
# Leaflet map
#m <- leaflet(width = "100%") %>%
#  addTiles() %>%
#  addPolygons(data = pdm, weight = 0.5, fillColor = "purple", strokeColor = "white", fillOpacity = 1, opacity = 1)
#m
m <- leaflet(width = "100%") %>%
  addTiles() %>%
  addPolygons(
    data = zoning, 
    stroke = TRUE, weight = 0.25, ,smoothFactor = 0.5, 
    color = "white", opacity = 1,
    fillColor = "grey",
    fillOpacity = 0.75,
    popup = zoning$zn_lclass)
m
st_write(zoning_int, dsn = "test.shp", driver = "ESRI Shapefile")
Writing layer `test' to data source `test.shp' using driver `ESRI Shapefile'
features:       297
fields:         5
geometry type:  Unknown (any)
GDAL Error 1: Attempt to write non-polygon (POINT) geometry to POLYGON type shapefile.
Failed to create feature 10 in test
Failed to create feature 10 in test
Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  : 
  Feature creation failed.
LS0tCnRpdGxlOiAiVG9wbyBNYXBwaW5nIFNwZWMgQ2hhbmdlIC0tIFRlc3RpbmciCm91dHB1dDogaHRtbF9ub3RlYm9vawpPcmdhbml6YXRpb246ICJDaXR5IG9mIFRvcm9udG8iCkRpdmlzaW9uOiAiSW5mb3JtYXRpb24gJiBUZWNobm9sb2d5IgpVbml0OiAiR2Vvc3BhdGlhbCBDb21wZXRlbmN5IENlbnRyZSIKQXV0aG9yOiAiUnlhbiBHYXJuZXR0IgpEYXRlOiAiTWFyY2ggMjQsIDIwMTgiCi0tLQoKIyMgMS4wIFByZXJlcXVlc2l0cwoKCmBgYHtyfQojIFJlYWRMb2FkIGxpYnJhcmllcyBhbmQgbG9hZCBzaGFwZWZpbGUgd2l0aCBzZiBwYWNrYWdlCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHNmKQpsaWJyYXJ5KHJhc3RlcikKbGlicmFyeShsZWFmbGV0KQoKc2V0d2QoIn4vR2l0aHViL2RhdGFTdG9yaWVzL3RvcG9TcGVjIikKCnpvbmluZyA8LSBzdF9yZWFkKCJ6b25pbmdfc3Vic2V0X3dnczg0LnNocCIpCnBkbSA8LSBzdF9yZWFkKCJwZG1fc3Vic2V0X3dnczg0LnNocCIpCmBgYAoKCmBgYHtyfQpwZG0gPC0gc2VsZWN0KHBkbSwgZ2VvbWV0cnkpCgpoZWFkKHBkbSkKYGBgCgoKYGBge3J9CiMgSW50ZXJzZWN0IHpvbmluZyB3aXRoIHBkbSBib3VuZGFyaWVzCgp6b25pbmdfaW50IDwtIHN0X2ludGVyc2VjdGlvbih6b25pbmcsIHBkbSkKCmNvdW50KHpvbmluZ19pbnQkem5fbGNsYXNzKQpoZWFkKHpvbmluZ19pbnQpCnBsb3Qoc3RfZ2VvbWV0cnkoem9uaW5nX2ludCksIGNvbCA9IHNmLmNvbG9ycyg0LCBjYXRlZ29yaWNhbCA9IFRSVUUpLCBib3JkZXIgPSAnZ3JleScsIGF4ZXMgPSBUUlVFKQoKCiNwbG90KHN0X2dlb21ldHJ5KG5jKSwgY29sID0gc2YuY29sb3JzKDEyLCBjYXRlZ29yaWNhbCA9IFRSVUUpLCBib3JkZXIgPSAnZ3JleScsIAojICAgIGF4ZXMgPSBUUlVFKQpgYGAKCgpgYGB7cn0KCmBgYAoKCmBgYHtyfQojIENvbnZlcnQgdGV4dCBhdHRyaWJ1dGVzIHRvIGludGVnZXIKCnpvbmluZ19pbnQgPC0gbXV0YXRlKHpvbmluZ19pbnQsIHpuX251bSA9IGlmZWxzZSh6bl9sY2xhc3MgPT0gIkNvcmUiLCAxLAogICAgICAgICAgICAgICAgICAgICAgaWZlbHNlKHpuX2xjbGFzcyA9PSAiSW5kdXN0cmlhbCIsIDIsCiAgICAgICAgICAgICAgICAgICAgICBpZmVsc2Uoem5fbGNsYXNzID09ICJSdXJhbCIsIDMsCiAgICAgICAgICAgICAgICAgICAgICBpZmVsc2Uoem5fbGNsYXNzID09ICJTdWJ1cmJhbiIsIDQsCiAgICAgICAgICAgICAgICAgICAgICAwKSkpKSkKYGBgCgoKYGBge3J9CgpwYWxldHRlIDwtIGNvbG9yQmluKGMoJyM3YjMyOTQnLAogICAgICAgICAgICAgICAgICAgICAgJyNjMmE1Y2YnLAogICAgICAgICAgICAgICAgICAgICAgJyNhNmRiYTAnLAogICAgICAgICAgICAgICAgICAgICAgJyMwMDg4MzcnKSwgCiAgICAgICAgICAgICAgICAgICAgIGJpbnMgPSBjKDEsIDIsIDMsIDQpKQpgYGAKCgpgYGB7cn0KIyBMZWFmbGV0IG1hcAoKI20gPC0gbGVhZmxldCh3aWR0aCA9ICIxMDAlIikgJT4lCiMgIGFkZFRpbGVzKCkgJT4lCiMgIGFkZFBvbHlnb25zKGRhdGEgPSBwZG0sIHdlaWdodCA9IDAuNSwgZmlsbENvbG9yID0gInB1cnBsZSIsIHN0cm9rZUNvbG9yID0gIndoaXRlIiwgZmlsbE9wYWNpdHkgPSAxLCBvcGFjaXR5ID0gMSkKI20KCgptIDwtIGxlYWZsZXQod2lkdGggPSAiMTAwJSIpICU+JQogIGFkZFRpbGVzKCkgJT4lCiAgYWRkUG9seWdvbnMoCiAgICBkYXRhID0gem9uaW5nLCAKICAgIHN0cm9rZSA9IFRSVUUsIHdlaWdodCA9IDAuMjUsICxzbW9vdGhGYWN0b3IgPSAwLjUsIAogICAgY29sb3IgPSAid2hpdGUiLCBvcGFjaXR5ID0gMSwKICAgIGZpbGxDb2xvciA9ICJncmV5IiwKICAgIGZpbGxPcGFjaXR5ID0gMC43NSwKICAgIHBvcHVwID0gem9uaW5nJHpuX2xjbGFzcykKbQoKYGBgCgoKYGBge3J9CnN0X3dyaXRlKHpvbmluZ19pbnQsIGRzbiA9ICJ0ZXN0LnNocCIsIGRyaXZlciA9ICJFU1JJIFNoYXBlZmlsZSIpCgpgYGAKCgo=